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SUMMARY 

The development of an internal layer in a turbulent boundary layer 
flow over a curved hill is investigated numerically. The turbulence field 
of the boundary layer flow over the curved hill is compared with that of a 
turbulent flow over a symmetric airfoil (which has the same geometry as the 
curved hill except that the leading and trailing edge plates were removed) 
to study the influence of the strongly curved surface on the turbulence 
field. The turbulent flow equations are solved by a control-volume based 
finite difference method. The turbulence is described by a 
multiple-time-scale turbulence model supplemented with a near-wall 
turbulence model. Computational results for the mean flow field (pressure 
distributions on the walls, wall shearing stresses and mean velocity 
profiles), the turbulence structure (Reynolds stress and turbulent kinetic 
energy profiles), and the integral parameters (displacement and momentum 
thicknesses) compared favorably with the measured data. Computational 
results show that the internal layer is a strong turbulence field which is 
developed beneath the external boundary layer and is located very close to 
the wall. Development of the internal layer was more obviously observed in 
the Reynolds stress profiles and in the turbulent kinetic energy profiles 
than in the mean velocity profiles. In this regard, the internal layers is 
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significantly different from wall-bounded simple shear layers in which the 
mean velocity profile characterizes the boundary layer most 
distinguishably . Development of such an internal layer, characterized by an 
intense turbulence field, is attributed to the enormous mean flow strain 
rate caused by the streamline curvature and the strong pressure gradient. 

In the turbulent flow over the curved hill, the internal layer begin to 
form near the forward corner of the hill, merges with the external boundary 
layer, and develops into a new fully turbulent boundary layer as the fluid 
flows in the downstream direction. For the flow over the symmetric airfoil, 
the boundary layer began to form from almost the same location as that of 
the curved hill, grew in its strength, and formed a fully turbulent 
boundary layer from mid-part of the airfoil and in the downstream region. 
Computational results also show that the detailed turbulence structure in 
the region very close to the wall of the curved hill is almost the same as 
that of the airfoil in most of the curved regions except near the leading 
edge. Thus the internal layer of the curved hill and the boundary layer of 
the airfoil were also almost the same. Development of the wall shearing 
stress and separation of the boundary layer at the rear end of the curved 
hill mostly depends on the internal layer and is only slightly influenced 
by the external boundary layer flow. 
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INTRODUCTION 


Turbulent flows subjected to various strain rates (in addition to the 
simple shearing strain rates) caused by streamline curvature, strong 
pressure gradient, separation and reattachment, swirl velocity and 
interaction of multiple number of shear layers are usually called "complex 
turbulent flows." The turbulence structure of such complex turbulent flows 
is more complicated than that of simple shear layer flows. Calculations of 
complex turbulent flows using various turbulence models such as the k-c 
turbulence models, algebraic Reynolds stress turbulence models (ARSM) , and 
Reynolds stress turbulence models (RSM) yield rather unsatisfactory 
computational results (Kline, Cantwell & Lilley 1982). Many turbulence 
models, improved by modifying the standard form turbulence equations 
(usually, the dissipation rate equation for k-e and ARSM and the 
pressure -strain correlation term for RSM) to yield better computational 
results for a few flow cases, have produced worse agreement with the 
measured data than the standard turbulence models for other classes of 
turbulent flows (Persen 1986). In recent years, a number of efforts have 
been made experimentally to better understand complex turbulent flows by 
separating the mechanisms that produce the extra strain rates. For example, 
experimental investigations of turbulent flows in curved channels (Gillis & 
Johnston 1983; So & Mellor 1973) and turbulent flows over a curved hill 
(Baskaran, Smits & Joubert 1987) were made to study the effect of 
streamline curvature on turbulence structure. 

Turbulent shear layers over curved surfaces are highly sensitive even 
to a small amount of streamline curvature (Bradshaw 1969 & 1973). Bradshaw 
(1969) proposed a curvature correction method based on an assumption that 
such turbulent flows can be characterized by a "curvature parameter", that 
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is, the ratio of boundary layer thickness (6) to radius of curvature (R) . 

In the curvature correction method, the mixing length is altered by a 
factor (F) given as; 

U/r 

F - 1 - 0 

JU/3n 

where U is the streamwise velocity, n is the coordinate normal to the 
streamline, r is the radius of curvature of the local streamline, and ft is 
a constant coefficient. Application of the method to a number of turbulent 
flows with varying streamline curvatures revealed that the constant 
coefficient (ft) needs to be adjusted for different flows to obtain 
computational results in good agreement with the measured data. To take 
into account the varying coefficient (0) for turbulent flows with different 
streamline curvature, a functional form of 0 which depends on the curvature 
Richardson number was also proposed (Gibson 1978). Many turbulence models 
incorporating a curvature correction method yield improved computational 
results for turbulent flows over mildly curved surfaces and the 
computational results help to better understand the turbulence structure of 
such flows; however, these turbulence models still fail to predict the 
turbulence field for turbulent flows with large streamline curvature. A 
discussion on the shortcomings of various curvature correction methods can 
be found in Gibson, Jones & Younis (1981). Gillis & Johnston (1983) found 
that turbulent flow over a curved surface reaches a "saturated state" as 
the streamline curvature is further increased and that many turbulence 
models did not yield meaningful results for such flows. Baskaran, Smits 6c 
Joubert (1987) presented experimental data for a turbulent boundary layer 
over a curved hill and a turbulent flow over a symmetric airfoil. The 
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curved hill and the airfoil shared the same geometry. They found that an 
internal boundary layer is formed beneath the external boundary layer over 
the curved hill and that the internal layer shares many similarities with 
the boundary layer flow over the symmetric airfoil. The experimental data 
also showed that the internal layer was insensitive to the external 
boundary layer and the curvature parameter (6/R). It was also suggested in 
Baskaran, Smits & Joubert (1987) that the effects of prolonged convex 
curvature can be more accurately represented by the ratio of the internal 
layer thickness to the surface radius of curvature than by the conventional 
curvature parameter. Here, detailed computational results for the turbulent 
boundary layer flow over the curved hill and the turbulent flow over the 
symmetric airfoil obtained using a multiple- time-scale turbulence model 
(Kim & Chen 1987; Kim 1988) are presented. The multiple- time- scale 
turbulence model (hereafter abbreviated as the M-S turbulence model), the 
near-wall turbulence model and the numerical method used herein are 
introduced below. 

It has been shown previously that the high Reynolds number M-S 
turbulence model yields accurate computational results for a number of 
complex turbulent flows. The wall- jet flow (Irwin 1973) and the 
wake-boundary layer interaction flow (Tsiolakis, Krause & Muller 1983) 
involve interaction of two turbulent shear layers. It is shown in Kim & 

Chen (1987) that the M-S turbulence model yields accurate computational 
results for these flows. It is also shown in Kim & Chen (1987) that the 
computational results for a backward- facing step flow (Kim, Kline & 

Johnston 1980) obtained using the M-S turbulence model compared somewhat 
more favorably with the measured data than those obtained using a Reynolds 
stress turbulence model. In the confined coaxial swirling jet (Roback & 
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Johnson 1983) , the axial velocity in the center region is retarded by the 
influence of the swirl velocity and a large reversed flow region is formed 
in the core region of the flow. Due to the existence of the reversed flow 
region, which is often referred to as a "flame stabilizer" in combustors, 
the flow is subjected to a large amount streamline curvature. This flow 
also contains a separated and reattaching shear layer in the corner region. 
It can be found in Kim & Chen (1987) that the M-S turbulence model yields 
significantly improved computational results compared with those obtained 
using the standard k-e turbulence models. The experimental study of a 
reattaching shear layer in a divergent channel (Driver & Seegmiller 1985) 
was primarily designed to study the effect of pressure gradient on the 
development of Reynolds stresses and the reattaching shear layer, to 
identify any deficiency in turbulence closure models, and thus to improve 
predictive capability of turbulence models. A number of turbulence models, 
such as the k-e turbulence models and algebraic Reynolds stress turbulence 
models (ARSM) , were shown to yield poor computational results for the flow 
(Driver & Seegmiller 1985) . It was also shown that a modified ARSM yielded 
computational results which are in good agreement with measured data. 
However, generality of the improved predictive capability of the turbulence 
model for other complex turbulent flows has not been shown yet. On the 
other hand, the computational results obtained using the M-S turbulence 
model compared somewhat more favorably with the measured data than those 
obtained using the modified ARSM. Recall that the turbulent transport of 
mass and momentum is governed by the time scale of the energy containing 
large eddies and the dissipation of the turbulent kinetic energy is 
governed by the time scale of the fine scale eddies (Lumley 1983). In M-S 
turbulence models (Hanjelic, Launder & Schiestel 1980; Kim & Chen 1987; Kim 
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1988), the turbulent transport of mass and momentum is described using the 
time scale of the large eddies and the dissipation rate is described using 
the time scale of the fine-scale eddies. In this regard, a well posed M-S 
turbulence model can describe the turbulent transport process better than 
the single- time-scale turbulence models such as the k-e, ARSM, and RSM 
turbulence models. Single- time-scale turbulence models yield reasonably 
accurate computational results for simple turbulent flows; however, the 
predictive capability degenerates rapidly as turbulent flows to be solved 
become more complex. This nature may due to the use of a single time scale 
to describe both the turbulence transport mechanism and the dissipation 
rate. 

In numerical calculations of turbulent flows, wall function methods 
are most frequently used to model the near-wall region. The wall function 
methods have been derived from a logarithmic velocity profile which usually 
prevails in wall-bounded simple shear layer flows. The wall function 
methods are not valid if the logarithmic velocity profile no longer 
prevails in the near-wall region. For the turbulent flows considered 
herein, a strong inequilibrium turbulence field is formed in the region 
very close to the wall. Such a turbulence field can not be described 
adequately by the wall function methods . Many other cases for which the 
wall function methods become invalid can be found in Kim (1988). A detailed 
discussion on various approaches to be used in place of wall function 
methods, including their advantages and disadvantages, can also be found in 
Kim (1988). In the near wall turbulence model used herein, the turbulent 
kinetic energy equations are extended to include the near-wall low 
turbulence region and the energy transfer rate and the dissipation rate 
inside the near-wall layer are obtained from algebraic equations. The 
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algebraic equations were obtained from a k- equation turbulence model 
(Wolfshtein 1969). It would be appropriate to classify the method as a 
"partially low Reynolds number approach" to distinguish it from other 
classes of methods. This approach was first used in Chen and Patel (1987) 
to solve turbulent flows over airfoils. Recall that the turbulence length 
scale in the region very close to a wall is related to the normal distance 
from the wall while that of external flows is related to the flow field 
characteristics (Roshko 1976). This characteristic of wall bounded 
turbulent flows can be described quite accurately by the "partially low 
Reynolds number turbulence models." Development of the near-wall turbulence 
model and its application to fully developed turbulent channel and pipe 
flows can be found in Kim (1988) . It has been shown in the reference that 
the present near- wall turbulence model can resolve the over- shoot phenomena 
of the turbulent kinetic energy and the dissipation rate in the region very 
close to the wall and that significantly improved computational results for 
the turbulence structure in the near-wall region are obtained. It is also 
interesting to note that a similar k-equation turbulence model, which forms 
the basis of the present near -wall turbulence model, yields fairly accurate 
computational results for a fully developed unsteady turbulent pipe flow 
(Tu & Ramaprian 1983). Incorporation of the same near-wall turbulence model 
into a k-c turbulence model and its application to a yet different class of 
complex turbulent flows such as a supersonic turbulent flow over a 
compression ramp and a transonic flow over an axisymmetric curved hill can 
be found in Kim (1989a; 1989b). 

The numerical method used herein is based on the pressure correction 
method (Patankar 1980) which has been used most extensively to solve 
incompressible flows. However, the present numerical method is applicable 
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for both incompressible and compressible flows with arbitrary, complex 
geometries. The capability to solve compressible flows is achieved by 
including a convective incremental pressure term into the pressure 
correction equation (Kim 1989a) . The accuracy and the convergence nature of 
the numerical method have been demonstrated by solving a number of flow 
cases. The example problems considered in Kim (1989a; 1989b) include: 
developing channel and pipe flows, a two-dimensional laminar flow in a 90 
degree bent channel, polar cavity flows, a turbulent supersonic flow over a 
compression ramp, and a shock wave - turbulent boundary layer interaction 
in transonic flow over a curved hill. 

TURBULENCE EQUATIONS 

As introduced in the previous section, the M-S turbulence model yields 
accurate computational results for a number of complex turbulent flows 
subjected to various different types of extra strains. However, the model 
is relatively new and is discussed in some detail below. 

The M-S turbulence model is based on a simplified split- spectrum 
method (Hanjelic, Launder & Schiestel 1980). The transport equations for 
the turbulent kinetic energy of large eddies (kp) and that of fine scale 
eddies (k t ) are given as; 
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where Pr is the production rate of turbulent kinetic energy, and Cp and e t 
are the energy transfer and dissipation rates, respectively. Eqs . (1-2) 
imply that the turbulent kinetic energy is generated by the instability of 
the mean fluid motion, is transferred to the high wave number region, and 
is dissipated by the molecular viscosity of the fluid. This mathematical 
model is consistent with the physically observed evolution of the turbulent 
kinetic energy (Tennekes & Lumley 1972) except that the cascade process of 
the turbulent kinetic energy is over-simplified and is represented by the 
single energy transfer rate. This over- simplification is still better 
justified than the single- time-scale turbulence models if one considers 
that only the generation and dissipation of turbulent kinetic energy are 
considered in the latter classes of turbulence models. 

The energy transfer rate and the dissipation rate away from the 
near-wall region are given as; 

d€p 3 d£ p Pr 2 Pr£ p e p 2 

uj { O + ) } - c pl — + c p2 c p3 (3) 

axj 3xj a ep 3xj kp kp kp 

3e t 3 v t de t e p 2 € p £ t € t^ 

uj - { O + ) } - c t i + c t 2 - c t3 (4) 

3xj 3xj a et 3y k t k t k t 

where the load functions in eqs. (3-4) were obtained from a physical 
dimensional analysis, and Cp^ and c t ^ (i«*l,3) are the model constants. 

These model constants were obtained by solving a five by five system of 
equations and by numerically optimizing one model constant (c t ^) to yield 
the best solution for a fully developed channel flow (Laufer 1949) and a 
plane jet exhausting into a moving stream (Bradbury 1965). One equation for 
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the model constants has been obtained from the equilibrium turbulent flow 
condition. Two equations have been obtained by transforming the 
multiple- time -scale turbulence equations into asymptotic turbulence growth 
rate equations which are equivalent to that of Harris, Graham & Corrsin 
(1977). The other two equations have been obtained by transforming the 
present turbulence equations into asymptotic turbulence decay rate 
equations which are equivalent to that given in Harlow and Nakayama (1968). 
The turbulence model constants are given as; ^^=0.75, a 6 p=1.15, u^ t “0.75, 
cr ct -1.15. Cp^-0.21, Cp2“1.24, Cp 3 - 1.84, c^-0.29, c t 2 “ 1.28, and c t3 «1.66. 
Further discussion on the establishment of these turbulence model constants 
can be found in Kim & Chen (1987). 

The eddy viscosity away from the near-wall region is given as; 

f k2 / e P - c /i k2/e t 

where c^f (*=0.09) is a constant coefficient and c^ - c /if c t/ e p t ^ ie 
effective eddy viscosity coefficient. In this equation, the turbulence 
length scale is related to the energy transfer rate of the energy 
containing large eddies rather than the dissipation rate of the fine-scale 
eddies. Experimental data show that the effective eddy viscosity 
coefficient (c^) decreases as the ratio Pr/e t increases and increases as 
the same ratio decreases. The present M-S turbulence model yields the 
solution in such a way that: Pr>€p>« t and c^<c^f for P r >€ t ; P r “€p-c t and 
c /i“ c /if f° r p r“ e t» and c p <c t and c /i >c /if por p r <<c t* Thus the eddy viscosity 
given as eq. (5) yields a variation of the c^ which is in good agreement 
with that observed in experimental data as well as that of the generalized 
algebraic stress turbulence models (Rodi 1972; Launder 1982; Kim & Chen 
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1988) . 


The eddy viscosity equation proposed by Hanjelic et al . , which is 
given as ^t* c /i(^ c p + ^ c t)^ c p/ € p » compatible with near wall mixing length 
theory or with the standard near -wall analysis only when k t vanishes in the 
near wall region. On the other hand, the eddy viscosity given in eq, (5) is 
compatible with both the near wall mixing length theory and the wall 
function method. 

The energy transfer rate and the dissipation rate inside the near-wall 
layer are given as; 



where e^c^fV^kV^/zcy the standard dissipation rate for near-wall 
turbulent flows in equilibrium state, f £ -l-exp( -A £ R t ) is a wall damping 
function, R t -k^/i/e^ is a turbulent Reynolds number, and A £ -c^£^/^/2 is a 
constant coefficient. The wall damping function was constructed in such a 
way that eq. (6) takes a limit value of 2u\a/y^ for y~0 and f € becomes unity 
in the region slightly away from the wall where the turbulence is in 
equilibrium state. For near-wall equilibrium turbulent flows, the 
production rate is approximately equal to the dissipation rate (e t ) and 
hence the energy transfer rate (cp) from the low wave number production 
range to the high wave number dissipation range has to be approximately 
equal to both of them. Recall that the production rate vanishes on the wall 
and grows to a peak value at the wall coordinate (y + ) approximately equal 
to 15. Hence eq. (6) may not be a good approximation for 0<y + <15. However, 
use of the vanishing boundary condition for the turbulent kinetic energy on 
the wall yields a growth rate of turbulent kinetic energy and a production 
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rate that are in good agreement with experimental data as well as with 
theoretical analysis (Kim 1988). 

The eddy viscosity in the near-wall layer is given as; 

"t " c /if ( 7 > 

e l 

where f^-l-l./exp(A^yR t + A2R^^) is a linear function of the distance from 
the wall in the viscous sublayer and becomes unity in the fully turbulent 
region. A^-0.025 and A2“0. 00001 have been used for the near wall layer (see 
Kim 1988). The eddy viscosity given as eq. (7) grows in proportion to the 
cubic power of the distance from the wall. It can be found in Kim (1988) 
that the near -wall analysis yields the same growth rate of the eddy 
viscosity in the region very close to the wall. The improved predictive 
capability of the present turbulence model is attributed to the use of the 
multiple- time-scales , the physically consistent eddy viscosity equation, 
and the near-wall turbulence model. 

NUMERICAL METHOD 

The incompressible turbulent flow equations are given as; 
d d 

— (pu) + — (pv) - 0. (8) 

dx dy 


d d d du d du dv dp 

— ( puu) + — ( puv ) « — (2p e — ) + — {/i e (— + — )) < 9 ) 

dx dy dx dx dy dy dx dx 
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d d d du 3v d 3v 3p 
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dx dy ax ay dx dy dy dy 


where x and y are spatial coordinates, u and v are cartesian velocities, p 
is the density, /i e is the effective viscosity, /i is the molecular 

viscosity, is the turbulent viscosity, and eqs. (8-10) follow from the 
conservation of mass, u-momentum, and v-momentum, respectively. In 
numerical calculation of incompressible flows, the conservation of mass 
equation is replaced by a pressure correction equation given as: 


d 

* a p'l 
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P K — 

ax 
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dy 
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where the superscript (*) denotes the current values, and k^ are 

coefficients relating the incremental velocities and the incremental 
pressure and the last term represents the mass imbalance. A similar 
pressure correction equation for both incompressible and compressible flows 
and its derivation can be found in Kim (1989a). In the present numerical 
method, the velocities are located at the same grid points and the pressure 
is located at the centroid of the cell formed by the four adjacent velocity 
grid points. Any other flow variables such as the turbulent kinetic energy 
and the dissipation rate are located at the velocity grid points. 

The control volume for the pressure correction equation is defined as 
the cell enclosed by the four neighboring grid points. The 
velocity-pressure decoupling is eliminated by treating the pressure 
correction equation as a continuous form partial differential equation 
(hereafter, abbreviated as p.d.e.) rather than treating it as a constraint 
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condition at each grid point. In the former case, the discrete pressure 
correction obtained from eq. (11) becomes a five-diagonal system of 
equations for rectangular grids. The off-diagonal terms are contributed 
only by the non- orthogonal grids and their magnitude is much smaller than 
that of the diagonal terms even for highly skewed grids. On the other hand, 
in the latter method (Vanka, Chen & Sha 1980), the discrete pressure 
correction equation, obtained by directly substituting the incremental 
pressure - incremental velocity relations into the conservation of mass 
equation, yields a nine-diagonal system of equations for orthogonal grids. 
This discrete pressure correction equation can yield a velocity-pressure 
decoupled solution, whereas the former equation does not. The present 
pressure correction equation is more like that of Braza, Chassaing & Ha 
Minh (1986) than the more standard pressure correction equation (Vanka, 

Chen & Sha 1980) in the sense that the equation is dealt as a continuous 
form p.d.e. In solving the discrete system of equations, the off-diagonal 
terms are moved to the load vector term and the resulting system of 
equations can be solved using a tri-diagonal matrix algorithm (TDMA) . 
Incorporating the off -diagonals into the load vector term did not 
degenerate the convergence rate (Kim 1989a) . 

In control -volume methods, the discrete system of equations is derived 
by integrating the governing differential equations over the control 
volume. For curvilinear grids, the required number of interpolations to 
obtain flow variables at the cell boundaries is significantly reduced by 
using the present grid layout. The enhanced convergence rate, even when 
highly skewed, unequally spaced, curved grids are used, is attributed to 
the pressure correction equation cast in a continuous form p.d.e. and the 
grid layout which required fewer interpolations. As shown in this paper, 
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the method yields accurate computational results (i.e., free of numerical 
wiggles) for grid aspect ratios as large as a few thousand. The pressure 
correction equation cast in a continuous form p.d.e. is undoubtedly more 
advantageous than the more standard pressure correction equation (Vanka, 
Chen & Sha 1980) for the reasons discussed above. 

COMPUTATIONAL RESULTS 

The turbulent flow over a curved hill is schematically shown in Figure 
1. The chord length (C) and the height (H) of the curved hill and the 
airfoil are 1.284 meters and 0.208 meters, respectively. The unit Reynolds 
number based on the free stream velocity (^-20 m/sec) is 1 . 33xl0^/meter . 

In external flows, the location of the far field boundary and the boundary 
conditions specified at the far field can influence the computational 
results. In the present study, the dependence of the computational results 
on the grid size, the location of the far field boundary, and the boundary 
conditions prescribed at the far field boundaries has also been studied. 

In the experimental investigation of the turbulent flow over the 
curved hill, the flow was made turbulent using a trip wire located 0.65C 
upstream of the forward corner of the curved hill. Hence in the numerical 
calculations, the upstream boundary (30^) was located 0.61C upstream of the 
forward corner of the curved hill. Inclusion of the trip wire in numerical 
calculation of the entire flow field is prohibitive at present due to the 
limitation imposed by the computational resources. Thus, the inlet boundary 
conditions for the tangential velocity, the turbulent kinetic energies, and 
the dissipation rates (Cp and c t ) were obtained from experimental data for 
a fully developed boundary layer flow over a flat plate (Klebanoff 1955). 
The non-dimensional velocity and the turbulent kinetic energy profiles were 


16 



scaled to yield a boundary layer thickness of 0.044 meters at the inlet 
boundary. This inlet boundary condition is somewhat different from that of 
the tripped turbulent flow; however, it is considered to be a reasonable 
approximation to the tripped turbulent flow since development of the 
internal layer is less dependent on the external boundary layer flow as can 
be found in the following discussion. In a series of numerical tests for 
the flow over the curved hill, the downstream boundary was located at 1.4C 
and 1.9C downstream of the rearward corner of the hill, and the top 
boundary was located at 9.3H, 12. 3H, and 15. 9H away from the wall, 
respectively. The computational results were not sensitive to the location 
of the downstream boundary as long as it is located reasonably far away 
from the reversed flow region. Location of the top boundary (d^) 
influenced the computational results substantially. For the latter two 
cases (e.g., 12. 3H and 15. 9H), the computational results in the near-wall 
region were virtually the same and those in the far region differed by no 
more than one percent. The computational results presented herein were 
obtained with dfy a nd 3^3 located at 15. 9H and 1.4C away from the wall and 
the rearward corner, respectively. The boundary conditions tested are: 

(1) Dirichlet boundary condition (DBC) for u and v along dfy'* and 
vanishing gradient boundary condition (NBC) for u and v and uniform 
pressure along 603 

and 

(2) NBC for u and v and uniform pressure along 3^2 and 303 . 

In each case, u=V“kp=k t =0 was specified on the solid wall boundary ( 3 Q 4 ) 
and the NBC for kp , k t , Cp, and e t was specified along 3ft2 and #^3 • The 
computational results obtained using the two different sets of boundary 
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conditions were virtually the same when the external boundary and 3 O 3 ) 

was located far enough so that the extent of the external boundary does not 
influence the computational result. A part of these computational results 
are also shown in the following. 

For the turbulent flow over the symmetric airfoil, the far field 
boundaries were initially separated by the same distance from the airfoil 
as in the curved hill flow. Then the distances of the upstream (30^) and 
the downstream (<^ 3 ) boundaries were increased by approximately 0.5C, 
respectively, until the computational results became independent of the 
extent of the external boundary. The computational results presented herein 
were obtained with the upstream (30^) and the downstream ( 3 Q 3 ) boundaries 
separated by 2 . 1C and 3 . 1C from the forward and the rearward corners of the 
airfoil, respectively . Uniform flow velocities were prescribed at the 
upstream boundary. In the experiment, the flow was made turbulent using a 
trip wire installed at the leading edge region of the airfoil. However, 
inclusion of the trip wire in the numerical calculation of the entire flow 
is not feasible for the reason discussed previously. Hence in the numerical 
calculation, a small amount of turbulent kinetic energy (k- 2 . SxlO’^U^) and 
a dissipation rate (c t ~40 m^/sec^) were prescribed at the inlet boundary so 
that the flow became fully turbulent from the leading edge region of the 
airfoil. For this boundary condition, the turbulent viscosity at the inlet 
boundary is approximately equal to the kinematic viscosity of the fluid. 

The boundary conditions along the 8^2 an< * ^3 boundaries were the same as 
the first case of the curved hill flow. A symmetric boundary condition was 
used for all flow variables along the symmetry line in front of and behind 
the airfoil. Boundary conditions along the airfoil were the same as those 
used in the flow over the curved hill. 
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A number of meshes were also tested in conjunction with the domain 
independence study to obtain grid independent solutions. In each case, the 
discrete finite difference system of equations was solved iteratively until 
the relative error for all flow variables became less than the prescribed 
convergence criterion (5X10~~*). A 143x77 mesh for the flow over the curved 
hill and a 146x77 mesh for the flow over the airfoil are shown in Figure 2. 
The grid spacings transverse to the flow direction are the same in each 
case. The average grid size (Ay + ) to the first grid point from the wall is 
approximately 0.8 (i.e., An-2xl0’^ meters). The grid size in the normal 
direction was increased by a factor of 1.15 approximately. The partition 
between the near-wall layer and the external region was located at y + «200 
and 27 grid points were used inside the near-wall layer. The M-S turbulence 
model predicted the reversed flow regions at the rear end of the curved 
hill and the airfoil for all the cases tested in the present study. All the 
computational results presented herein, except the streamline contours 
shown in Figure 3, were obtained using the fine mesh. 

The streamline contours for the flows over the curved hill and the 
airfoil obtained using the medium mesh (92x60 grid points in each case) are 
shown in Figure 3. For the medium mesh, the average grid size (Ay + ) to the 
first grid point from the wall was 80 (i.e., An-0.001 meters). The 
partition between the near -wall layer and the external domain was located 
at y + ~300 and 3 grid points were allocated inside the near wall layer. The 
calculated mean velocity, pressure, and turbulence fields were in fair 
agreement with the measured data. 

The calculated static pressure distributions on the walls of the 
curved hill and the symmetric airfoil are compared with the measured data 
in Figure 4, where the pressure coefficients (Cp) were obtained by 
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normalizing the pressure distributions by the free stream dynamic pressure 
(O.SpU^). The calculated static pressure distribution along the symmetry 
line of the entire computational domain is also shown in the upper right 
corner of the figure. It can be seen from this figure that the numerical 
method does not yield an unphysical oscillatory solution for the mesh with 
the grid aspect ratio as large as a few thousand. Note that the measured 
data exhibit pressure jumps at the plateau of the hill and the airfoil. The 
pressure gradients caused by these jumps are believed to be strong enough 
to disturb the mean flow field and the disturbed mean velocity may leave 
some evidence in the measured wall shearing stress, at least slightly; 
however, such an evidence is not found in the measured wall shearing 
stress. It is not clear if these pressure jumps really existed or were 
caused by the interference of the wind tunnel wall or the measuring device. 
The cause of such pressure jumps was not clarified in Baskaran, Smits & 
Joubert (1987) . Such pressure jumps were not detected in the present 
calculations. It can be seen in the figure that the calculated pressure 
distributions compare favorably with the measured data, in general. The 
pressure distributions on the curved hill obtained using the two different 
boundary conditions are almost the same and almost collapse into a single 
line as shown in this figure. This indicates that the far field boundary 
condition can be implemented in one way or another as long as it is located 
far away from the wall. Also note that the measured pressure distribution 
on the wall of the airfoil is somewhat higher than the calculated result. 
This discrepancy might have been caused by the interference of the wall 
and/or measuring device in the experiments. This fact is further clarified 
in the following by comparing the measured and the calculated mean velocity 
profiles . 
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The calculated wall shearing stresses for the curved hill and the 
airfoil are shown in Figures 5 and 6, respectively, where the friction 
coefficients (cf) were obtained by normalizing the wall shearing stresses 
by the free stream dynamic pressure. As discussed in the previous 
paragraph, the wall shearing stresses obtained using the two different 
boundary conditions almost collapsed into a single line so that the present 
computational results are independent of the extent of the far field 
boundary and the boundary conditions prescribed at the far field boundary. 
The measured data and the computational results obtained using a curvature 
correction method (Baskaran, Smits & Joubert 1987) are also shown in this 
figure for comparison. It can be seen in Figure 5 that the calculated wall 
shearing stress for the curved hill is slightly smaller than the measured 
data near the inlet boundary. This discrepancy is attributed to the inlet 
boundary condition obtained from a fully developed turbulent boundary layer 
flow. Apparently, the mass flow rate entering the computational domain 
prescribed by the velocity profile for a fully developed turbulent flow is 
smaller than that for the developing velocity profile made turbulent with a 
trip wire. The influence of the slightly smaller mass flow rate is observed 
in the wall shearing stress distribution throughout the curved hill. In 
Figures 5 and 6, "S" represents the separation location and "R" represent 
the reattachment location. The M-S turbulence model predicted the small 
reversed flow region near the rear end of the curved hill. On the other 
hand, the calculated wall shearing stress using the curvature correction 
method is in close agreement with the measured data. This accurate 
computational result may due to the use of a boundary layer flow solver 
which incorporates the measured pressure distribution on the wall (or, 
equivalently, the external mean flow velocity). However, the curvature 
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correction method failed to predict the reversed flow region in the rear 
end of the curved hill. For the flow over the symmetric airfoil, the 
present computational result is in excellent agreement with the measured 
data. This excellent agreement is attributed to the the fact that the mass 
flow rate of the flat inlet velocity profile is in closer agreement with 
that of the experiment and to the use of the M-S turbulence model which can 
describe the turbulence field in a developing boundary layer flow more 
accurately. The wall shearing stress obtained using the same boundary layer 
flow solver and the curvature correction method does not compare very well 
with the measured data. This disagreement may be due to the curvature 
correction method which can not resolve the strong, developing turbulence 
field even though the measured pressure (or, the external potential 
velocity) is used in the boundary layer calculations. Again, the present 
computational result shows the reversed flow region located near the 
trailing edge of the airfoil while the curvature correction method fails to 
predict such a reversed flow region. It is shown in Figure 6 that the 
separation location for the airfoil compares more favorably with the 
measured data than it does in the case of flow over the curved hill. 

The mean velocity profiles for the flow over the curved hill and the 
airfoil at four downstream locations are compared with the measured data in 
Figure 7, where the measured mean velocity profiles in physical units were 
recovered from the measured velocity profiles given in non-dimensional form 
using the measured wall pressure. In this figure, s'-s-0.074 meters and 
0.074 meters represent the correction factor for effective origins of the 
two flows based on the pressure gradient parameter, see Baskaran, Smits & 
Joubert (1987) for details on the effective origin. It can be seen in this 
figure that the calculated and the measured mean velocity profiles exhibit 
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good comparison in general. For the flow over the airfoil, a fluid particle 
along the symmetry line is brought to zero velocity at the forward 
stagnation point and a strong adverse pressure gradient is formed in this 
region. Due to the strong adverse pressure gradient, the velocity in this 
region is significantly retarded. The computational results show that the 
mean tangential velocity for the airfoil in this region is smaller than 
that for the curved hill. The computational results also show that the 
tangential velocity for the airfoil is slightly larger than that for the 
curved hill in farther downstream locations. These results simply reflect 
that the inlet mass flow rate for the airfoil is larger than that for the 
curved hill in numerical calculations. Note that the measured mean velocity 
profiles for the airfoil at s' =1.139 meters and 1.862 meters are smaller 
than those for the curved hill. The smaller mean velocity for the airfoil 
in this region is due to the slightly higher pressure distribution in the 
plateau region of the airfoil which is subjected to a slight experimental 
uncertainty as discussed previously. 

The calculated displacement and momentum thicknesses are compared with 
the measured data in Figures 8 and 9, respectively. The displacement 
thickness ($*) and the momentum thickness ( 6 ) were calculated using 
equations given as : 
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where u p-Up W e^ n is the local potential velocity, Up W "“Uco(l' c pw^^^ is t ^ ie 
potential velocity on the wall, k is the curvature of the surface, and n is 
the normal distance from the wall. Eqs . (12) and (13) were used in So & 
Mellor (1973) and Baskaran, Smits & Joubert (1987) to study thin shear 
layers over curved surfaces, since the classical boundary layer theory and 
the relationship between the potential velocity (U^) and the wall pressure 
given as "p^+0 . SpU^-constant" does not hold very well for flows over 
curved surfaces. It can be seen in these figures that the calculated 
results and the measured data compare favorably in general. The calculated 
displacement and momentum thicknesses of the curved hill near the inlet 
boundary are slightly larger than the measured data. These slight 
discrepancies are again due to the inlet boundary condition as discussed 
previously. However, as the flow approaches the curved hill, these 
discrepancies disappear and the calculated results are in excellent 
agreement with the measured data. This observation partly shows that the 
development of the internal layer on the curved hill is only slightly 
influenced by the approaching external boundary layer flow. It has been 
shown in Figure 7 that the mean velocity profiles for the curved hill and 
the airfoil at s' *1.139 meters are almost identical. However, the 
displacement thickness and the momentum thickness of the flows over the 
curved hill and the airfoil, which were obtained from almost the same 
velocity profiles, are significantly different in the same region near the 
leading edges as shown in Figures 8 and 9. These disagreements may due to 
the use of the boundary layer theory for flows over curved surfaces 
somewhat beyond its applicable limit. Recall that eqs. (12) and (13) were 
originally intended to be used for flows over curved surfaces with very 
small curvature (So & Mellor 1973). The turbulent flow approaching the 
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curved hill is highly retarded due to the strong adverse pressure gradient 
existing near the leading edge of the curved hill and thus the displacement 
and the momentum thicknesses are increased significantly in this region. 

The same flow slightly beyond the leading edge is subjected to far stronger 
favorable pressure gradient and is accelerated enormously. Thus the 
displacement and the momentum thicknesses decrease abruptly. In farther 
downstream, the internal layer is formed gradually and thus these 
thicknesses grow gradually until the flow is subjected to separation at the 
rear end of the curved hill. Near the separation point, these integral 
parameters increase abruptly again. It was suggested in Baskaran, Smits & 
Joubert (1987) that the wavy nature of the displacement and momentum 
thicknesses of the flow over the curved hill might have been caused by the 
use of the trip wire. However, as discussed in the above, the present 
computational results suggest that the wavy nature of these integral 
parameters is inherent to the flow over the curved hill. It is also 
interesting to note that any turbulence model incorporating wall function 
methods may not be able to describe the turbulence field over the curved 
hill adequately partly due to the complex turbulence structure of the 
internal layer and partly due to the wavy nature of the boundary layer 
thickness. For example, the optimal distance from the wall where a wall 
function method can be applied is obscured due to the rapidly varying 
boundary layer thickness. 

The calculated Reynolds stress profiles for the flow over the hill and 
the flow over the airfoil at a number of downstream locations are shown in 
Figure 10. The same Reynolds stress profiles in wall coordinates are shown 
in Figure 11. It can be seen in Figure 10 that the calculated Reynolds 
stress profile at s=0.710 meters is slightly more spread out than the 
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measured data. This discrepancy is again attributed to the inlet boundary 
condition obtained from a fully developed boundary layer flow. At farther 
downstream locations, the calculated and the measures Reynolds stress 
profiles are in fair comparison with the measured data qualitatively and 
quantitatively. The form of the Reynolds stress profile at s-0.710 meters 
is similar to that of a wall-bounded simple shear layer flow and it belongs 
to the external boundary layer flow. It is shown in Figure 10 that the 
strength of the Reynolds stress of the external boundary layer flow decays 
gradually and that of the newly forming internal layer grows rapidly as the 
fluid travels in the downstream direction. At farther downstream locations, 
these two Reynolds stress profiles merge together and form a new profile 
which is similar to that of a wall-bounded simple shear layer flow. Thus 
the shape of the Reynolds stress profiles in most parts of the curved hill 
are distinctively different from those found in wall-bounded simple shear 
layer flows. It can be seen in Figure 11 that there exists a slight 
difference between the near-wall Reynolds stress profiles of the flow over 
the curved hill and the flow over the airfoil due to the different upstream 
turbulence structure. It can also be found that the maximum magnitude of 
the Reynolds stress is located at y + «20. The average physical distance from 
the wall to this location is approximately equal to 4x10'^ meters and thus 
detailed measurement of the turbulence structure in this region can be very 
difficult. Development of the maximum Reynolds stress along the flow 
direction is shown in Figure 12. For the flow over the curved hill, a large 
amount of Reynolds stress is convected downward from the upstream, while 
the magnitude of the Reynolds stress convected downward in the flow over 
the airfoil is almost negligible. However, as the flows approach the 
plateaus of the curved hill and the airfoil, the maximum Reynolds stresses 
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of these distinctively different flows become almost the same. Thus it is 
found that the highly intense turbulence field of the flow over the curved 
hill is almost independent of the upstream turbulence level and mostly 
dependent upon the curvature of the curved hill. 

The turbulent kinetic energy profiles of the curved hill and the 
airfoil at a number of downstream locations are shown in Figure 13 and 
those in wall coordinates are shown in Figure 14. Development of the 
maximum turbulent kinetic energy along the flow direction is shown in 
Figure 15. The turbulent kinetic energy profile at s-0.710 meters shows 
that the turbulence intensity of the tripped flow is larger than that of 
the fully developed boundary layer flow. However, at downstream locations, 
the calculated turbulent kinetic energy profiles are in good agreement with 
the measured data. This good agreement at downstream locations indicates 
that the development of the turbulence field in turbulent flows over 
strongly curved surfaces depends only slightly on the upstream turbulence 
intensities. It also indicates that the M-S turbulence model can correctly 
predict the turbulence fields in a strongly inequilibrium state. It is also 
shown in Figure 15 that the upstream turbulence intensity of the flow over 
the curved hill is higher than that of the flow over the curved hill near 
the inlet boundaries and that the turbulence intensities of these two flows 
become almost the same at farther downstream locations for the same reason 
as discussed previously. 


CONCLUSIONS AND DISCUSSION 

Detailed computational results for the turbulent flows over a curved 
hill and a symmetric airfoil obtained using a multiple- time* scale 
turbulence model were presented. The computational results show that an 
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internal layer is developed beneath the external boundary layer. The 
development of the internal layer depends mostly on the mean flow strain 
rates distributed over the curved hill and is only slightly influenced by 
the external boundary layer flow. The curved hill and the airfoil shared 
the same geometry. The pressure distributions on the walls, the wall 
shearing stresses, and the mean velocity profiles for the flow over the 
curved hill and the airfoil were almost the same. Thus the mean flow strain 
rates are also almost the same in these two flows. The internal layer of 
the flow over the curved hill, which was caused by almost the same mean 
flow strain rate as that of the airfoil, was also very similar to the 
boundary layer flow developed over the symmetric airfoil. It can be seen 
from these numerical results, as well as from the measured data, that the 
development of the turbulence field on strongly curved surfaces depends 
only slightly on the external (or, the upstream) flows and depends mostly 
on the mean flow strain rate caused by the curved surfaces. Thus the 
development of the turbulence fields in these two flows are similar to the 
saturated behavior observed in turbulent flows over strongly curved 
surfaces . 

In the turbulent flow over the curved hill, the pressure gradient 
changes from adverse to favorable and back to adverse again. The 
displacement and momentum thicknesses of the boundary layer flow subjected 
to this pressure gradient increased abruptly near the forward corner of 
the curved hill, decreased significantly and increased gradually along the 
curved hill, and then increased abruptly again near the rear end of the 
curved hill where the flow separates. Thus the wavy nature of these 
integral parameters as observed in the numerical results and in the 
measured data is considered to be the inherent nature of the flow over the 
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curved hill. 

The development of the internal layer was more obvious in the 
turbulence field (the Reynolds stress profiles and the turbulent kinetic 
energy profiles) than in the mean velocity profiles. Thus the internal 
layer can be characterized as a strong turbulence field developing 
underneath the external boundary layer and in the region very close to the 
wall. Development of such a strong turbulence field was caused by the mean 
flow strain rate. Therefore, development of the internal layer slightly 
lags that of the mean velocity field. 

In the flow over the curved hill, the turbulence intensity of the 
external boundary layer becomes weaker while that of the internal layer 
becomes stronger as the fluid flows in the downstream direction. Both the 
calculated and the measured data showed that the turbulence fields of the 
external boundary layer and the internal layer merged together in the 
middle region of the curved hill and formed a new turbulence field which is 
similar to that of a wall -bounded simple shear layer at the rear end of the 
curved hill. The internal layer was located very close to the wall and thus 
development of the wall shearing stress and the flow separation at the rear 
end of the curved hill were determined solely by the internal layer. The 
velocity profiles and the wall shearing stress distributions on the curved 
hill and the airfoil were almost the same. Thus the logarithmic velocity 
profiles, obtained by non-dimensionalizing the mean velocity profiles using 
the wall shearing stress, for the curved hill and the airfoil are also very 
similar as can be found in Baskaran, Smits & Joubert (1987). However, the 
Reynolds stress and the turbulent kinetic energy profiles slightly away 
from the near-wall region of the curved hill are significantly different 
from those of the airfoil. From these results, it is found that the 
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logarithmic velocity profile alone is not sufficient to characterize the 
complex turbulent flow over the curved hill. 
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FIGURE 1. - NOflENCLATURE FOR TURBULENT FLOW OVER A CURVED HILL. 


34 




(b) SYMMETRIC AIRFOIL. 


FIGURE 2. - DISCRETIZED FLOW DOMAINS. 
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FIGURE A. - STATIC PRESSURE DISTRIBUTIONS ON THE HALLS. FIGURE 5. - HALL SHEARING STRESS FOR CURVED HILL. 




FIGURE 6 . - HALL SHEARING STRESS FOR SYMMETRIC AIRFOIL. 


FIGURE 7. - MEAN VELOCITY PROFILES. 
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FIGURE 8. - DISPLACEMENT THICKNESS. 





FIGURE 11. - REYNOLDS STRESS IN THE NEAR-HALL REGION. 
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FIGURE 12. - MAXIMUM REYNOLDS STRESS DISTRIBUTIONS ALONG 
THE FLOW DIRECTION. 



0 0 0 0 0 0 .01 


FIGURE 14. - TURBULENT KINETIC ENERGY PROFILES IN THE 
NEAR-WALL REGION. 
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FIGURE 13. - TURBULENT KINETIC ENERGY PROFILES. 
NOTATIONS AS IN FIG. 10. 
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FIGURE 15. - MAXIMUM TURBULENT KINETIC ENERGY DISTRIBUTION. 
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